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Abstract. We present a novel numerical method, called Jacobi-predictor-corrector approach, 
for the numerical solution of fractional ordinary differential equations based on the polynomial inter- 
polation and the Gauss-Lobatto quadrature w.r.t. the Jacobi-weig ht function u(s) = (l-s) Q_1 (l+s) . 
j /' j ■ This method has the computational cost O(N) and the convergent order IN, where N and IN are, 

respectively, the total computational steps and the number of used interpolating points. The detailed 
error analysis is performed, and the extensive numerical experiments confirm the theoretical results 



and show the robustness of this method. 
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1. Introduction 

. This paper further discusses the numerical algorithm for the following initial value problem 

(N ■ 

D?x(t)=f{t,x(t)), a:< fc >(0) = 4*°. fc = 0,l,-- - ,\a] -1, (1.1) 
- ■ — i ■ 

| where a 6 (0, oo), can be any real numbers and D™ denotes the fractional derivative in the Caputo sense [15], 
■ defined by 



7S=5J /o(* " T) n -<*- l xW{T)dT, n - 1< a < n; 



I -^t 1 , a = n; 

where n := \a] is just the value a rounded up to the nearest integer, D™ is the classical nth-order derivative, 
and 

1 r* 

oD'^xit) = —— / (t - tY^x^cLt 
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is the Riemann-Liouville integral operator of order fi > 0. It is well known that the initial value problem (1.1) 
is equivalent to the Volterra integral equation [3,8,9] 



in the sense that if a continuous function solves (1.2) if and only if it solves (1.1). 

Many approaches have been proposed to reslove (1.1) or (1.2) numerically, such as [4,5,8-10]. Diethelm, 
Ford and their coauthors successfully present the numerical approximation of (1.2) using Adams-type predictor- 
corrector approach and give the corresponding detailed error analysis in [8] and [9] , respectively. The convergent 
order of their approach is proved to be min{2, 1 + a}, and the arithmetic complexity of their algorithm with steps 
N is 0(N 2 ), whereas a comparable algorithm for a classical initial value problem only give rise to O(N). The 
challenge of the computational complexity is essentially because fractional derivatives are non-local operators. 
This method has been modified in [4], where the convergent order is improved to be min{2, 1 + 2a} and 
almost half of the computational cost is reduced, but the complexity is still 0(N 2 ). There are already two 
typical ways which are suggested to overcome this challenge. One seems to be the fixed memory principle of 
Podlubny [15]. However, it is shown that the fixed memory principle is not suitable for Caputo derivative, 
because we cannot reduce the computational cost significantly for preserving the convergent order [8,10]. The 
other more promising idea seems to be the nested memory concept of Ford and Simpson [5,10] which can lead 
to 0(N log N) complexity, but still retain the order of convergence. However, the convergent order there cannot 
exceed 2. For the effectiveness of the short memory principle, in [10], a has to belong to the interval (0, 1); and 
in [5], a must be within (0,2). 

In this work, we apprehend the Riemann-Liouville integral from the viewpoint of a normal integral with a spe- 
cial weight function. Thus we can deal with it based on the theories of the classical numerical integration and of 
polynomial interpolation [16]. Then by using a predictor-corrector method, called Jacobi-predictor-corrector 
approach, we obtain a good numerical approximation to (1.2) with the convergent order IN, which is the num- 
ber of used interpolating points. Moreover, the computational complexity is reduced to O(N), the same as 
classical initial value problem, which is one of the most exciting and significant advantages of this algorithm. 

The organization of this paper is as follows. In Section 2, we present the Jacobi-predictor-corrector approach 
and its detailed algorithm. In Section 3, the error analysis of the numerical scheme is discussed in detail. The 
algorithm is simply modified in Section 4 to deal with the extreme cases. Two numerical examples are given 
in Section 5 to confirm the theoretical results and to demonstrate the robustness the algorithm. Concluding 
remarks are given in Section 6. Finally, some Jacobi-Gauss-Lobatto nodes and weights w.r.t. the Jacobi-weight 
functions (1 — s) a (1 + s) used in the numerical experiments are listed in Appendix at the end of this paper. 



In this section we shall derive the fundamental algorithm for numerically solving the initial value problems 
with Caputo derivative. It is the following transformation other than (1.2) itself that underlies this algorithm: 




(1.2) 



2. Jacobi-predictor-corrector approach 





(2.1) 



where 
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x(s) = x(^(l + s)), -1 < s < 1. 

We assume the function / to be such that a unique solution of (1.2) exists in some interval [0, T], specifically 
these conditions are (a) the continuity of / with respect to both its arguments and (6) a Lipschitz condition 
with respect to the second argument [7]. Thus by the theory of the classical numerical integration [16], we 
can approximate the integral in (2.1) using the Jacobi-Gauss-Lobatto quadrature w.r.t. the weight function 
uj(s) = (1 - s)"" 1 ^ + s)°. That is, 

k=0 ' ^ ' j=0 

where JN + 1, {ujj}^ , {sj}j= in (2.2) correspond to the number of, the weights of, and the value of the 
Jacobi-Gauss-Lobatto nodes in the reference interval [—1, 1], respectively. 

Let us define a grid in the interval [0, T] with N + 1 equi-spaced nodes U, given by 

ti =ih, i = Q,--- ,N, (2.3) 

where h — T/N is the step-length. Suppose that we have got the numerical values of x(t) at to, t\, • • • , t n , which 
are denoted as xq, xi, ■ ■ ■ , x n , separately (xq = x °^), now we are going to compute the value of x(t) at t n +i, 
i.e. x n+ i. 
By (2.2), 

M_1 t k It JN 

x(t n+ i)*i Y -^ x o k) + f^(-^) a Y u jfn+i(sj,x n+1 ( Sj )), (2.4) 

k=0 ' ^ ' j=0 

where 

/ n+1 ( s ,x„ +1 ( S ))=/(^i(l+ S ),x(^i(l + S ))), -1< S <1; 

x n+1 (s) = x(^(l + «)), -1 < s < 1. 

To do the summation of the second term of (2.4), first we need to evaluate the value of / at the point t n+ i 
since fn+i(s.JN,x n +i(sjN)) = f(t n +i, x{t n +i)) ■ This value can be numerically got by using the interpolation 
of / w.r.t. t based on the known values of / at the equi-spaced nodes which are in the "neighborhood" of t n+ \. 
For the other values of f n +\ (sj, x n +i(sj)) , < j < JN — 1, we can also obtain them based on the interpolation 
of / at the equi-spaced nodes located in the "neighborhood" of Sj (should be (1 + s 3 )t n+ i/2 as to variable t). 
Denote IN as the number of equi-spaced nodes used for the interpolation. To start this procedure, the values 
XqiXx,- ■■ ,xin-i should be known first, and should be accurate enough for not deteriorating the accuracy of 
the algorithm. 

Now we make it clear what the "neighborhood" means. For predicting the value of / at t n+ i, we use the 
values of / at t n -iN+ii tn-iN+2, • • ■ , tn-i> t n - For getting the values of f n +i at Sj, < j < JN — 1, the way 
to choose the "neighborhood" equi-spaced nodes is as follows: (i) try to make \ln — rn\ as small as possible; (ii) 
make In > rn if possible, where In and rn are respectively the number of equi-spaced points used on the left 
and right hand side of Sj (should be (1 + Sj)t n+ i/2 as to variable t), obviously, In + rn = IN. So far, we arrive 
at the predictor-corrector formulas of (2.1) as 

W-l t k ^ 1 t n ! a ( JN ~ 1 

x n +i= J kT x ° + Wo) ( J 2^) Q ( E "jfn+lj +U) JN f(t n+ i,a£ +1 )), (2.5) 

k=0 ' ^ > j=0 
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and 

= E ^ + fL(Y) a S w ^' (2 - 6) 

where {fn+i j}j=o m (2-6) means that all the values of f n +i at the Jacobi-Gauss-Lobatto nodes are got by using 
the interpolations based on the values of {f(ti, Xi)}"_ ; whereas {/n+i,j}/=o m (^-5) are obtained by using 
the interpolations based on the values of {f(ti,Xi)}f_ and f(t n +i,x^ +1 ). The algorithm for realizing (2.5) and 
(2.6) is detailedly described as: 

Step 0. Some notations: 

sum := saving the value of the second summation in the right hand of (2.6) (or (2.5)), 
initially sum = 0; 

Pl{t) := the interpolating polynomial passing through (to, f(to, xq)), {tit f{tii x i))i ' ' ' i 

(tlN-l, f(tlN-l,XlN-l))', 

Pr n +i{t) := the interpolating polynomial passing through (t n - IN+ i, f(t n - IN+ i, x n -iN+i)), 

(tn-IN+2, f{t n -IN+2, X n -IJf+<l)), ■ ■ ■ , (t n , f(t n , X n ))] 

Cr n+ i(t) := the interpolating polynomial passing through (t n - IN+2 , f(t n - IN+2 , x n -i N+2 )), 

(tn-IN+3, J '(tn-IN+3, X n - IN+3 )), ■ • • , (t n+ i, f(t n+ i, X^ +1 ))\ 

In := \IN/2], to evaluate f n+ i(sj, i n +i(sj)J, the expected number of the interpolating 

equi-spaced nodes on the left hand side of Sj ; 
rn :— |_J"iV/2j , to evaluate f n +i(sj, x n +i(sj)^) , the expected number of the interpolating 

equi-spaced nodes on the right hand side of Sj ; 
le := the number of the interpolating equi-spaced nodes that can be used on the left of Sj] 
P n +i{t) := the interpolating polynomial passing through (t{ e _{ n > f(tle-ln, Xle-ln))> 

— %le— ln-\-l 

Step 1. To start the procedure: 

Compute Xi,X2, • ■ ■ ,xin-i by a single step method (e.g., the Improved- Adams' methods in [4]) with a suf- 
ficiently small step- length ho such that Xi, i — 1, 2, • • • , IN — 1, are accurate enough for not deteriorating the 
accuracy of the method we are discussing. 

Step 2. To predict: 

do j = 0, • • • , JN 

if le < In (the number of the equi-spaced nodes located on the left hand of Sj (should be (1 + Sj)t n+ i/2 as 
to variable t) is equal to / less than what we expect) 
sum — sum + ujjPl((l + Sj)t n +i/2J 

else if le + rn > n + 1 (the number of the equi-spaced nodes located on the right hand of Sj (should be 
(1 + sj)t n+ i/2 as to variable t) is equal to / less than what we expect) 
sum — sum + WjPr n+ i((l + Sj)t n +i/2) 
else 

sum — sum + uijP n+ i ((1 + Sj)t n+ i/2) 
enddo 

p _ v r«]-l + i (fc) 1 (t n + 1 \a 

x n+l — 2^k=0 fc! x + T(a) V 2 J iUm ' 
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Step 3. To correct: 

do j = 0, • • • , JN - 1 
if le < In 

sum = sum + uijPlUl + Sj)i n+1 /2) 

else if le + rn > n + 2 

sum = sum + UjCr n+ i ((1 + Sj)t n +x/2\ 

else 

sum — sum + uijP n+ i ((1 + Sj)t n +i/2) 
enddo 

x n +i = EIIo" 1 HhT^ + T^)( h f 1 ) a ' {sum + uj JN f{t n+1 ,x^ +1 )). 
We call the above algorithm J acobi-predictor- corrector approach. 

Although the description of this algorithm seems tedious, its operation is simple and mechanical. It can be 
observed that for the computation of x n +i, only changeless 2(JN + 1) values are needed, each of which should 
be interpolated by changeless IN nearby values; whereas in [4,8], it should take 0(n + 1) multiplications and 
divisions. In other words, the computational cost here has no relation with the variable n + 1, just depends 
on the number of the interpolating nodes IN and the number of Jacobi-Gauss-Lobatto nodes JN + 1, so, 
to approximate x(T), the total computational cost is O(N), comparing with 0(N 2 ) in [4,8] and O(NlogN) 
in [5, 10], which is one of the most significant advantages of this algorithm. 

3. Error Analysis 

First, we introduce four notations that will be used in the following error analysis. The piecewise interpolating 
polynomial based on the IN nodes of { (t i; f(ti, Xi)) } is denoted by AI IN ' P / n+ i(s), where — 1 < s < 1; the 
one based on the IN nodes of { (ti, ffa, x^) }™ =Q and f{t n +i, is wrote as AI IN f n+ i(s), where 

— 1 < s < 1; the one based on the IN nodes of { (ti, f{ti, x(U))) }"_ Q is signified by PI IN P f n+1 [s), where 

— 1 < s < 1; the one based on the IN nodes of { (ti, f(ti, x(ti))} }™^ Q is denoted as PI IN f n+ i(s), where 
— 1 < s < 1. Note that Xi is the numerical solution and x(U) is the exact solution. 

Here authors state that the idea of the error analysis below is inspired by that in Kai Diethelm, Neville J. 
Ford and Alan D. Freed's paper [9]. 

3.1. Some preliminaries and a useful lemma 

Let <jj a 'P(x) = (1 — x) a (l + x)P , a >— 1,/3>— lbea Jacobi- weight function in the usual sense. As illustrated 
in [1,11-14,17,18], the set of Jacobi polynomials { J™ ,/3 (x)}^ =0 forms a complete (—1, l)-orthogonal system, 
where L 2 ^^ (—1, 1) is a weighted space defined by 

1,1) = {v : v is measurable and || v \\ ua ,fi< oo}, (3.1) 

equipped with the norm 




and the inner product 

(u,v) u o,,3= J u(x)v(x)uj a ^ (x)dx. (3-3) 

For bounding the approximation error of Jacobi polynomials, we need the following non-uniformly-weighted 
Sobolev spaces as in [17]: 

iC..„,,(-l, 1) : = : 9 k x v S Ll a+hifi+k (-l, 1), < k < m}, 
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equipped with the inner product and the norm as 

m 

{u,v) m>u c,f> t * = ^2(d*u,d%v) ua +k,i3+k, (3.4) 



k=0 

and 



II V \\m,u,".P,*= y / (« J «)m^°./3,*- (3-5) 

For any continuous functions u and v on [—1,1], we dehne a discrete inner product as 

N 

(u,v) N =^2u(x j )v(x j )uj j , (3.6) 

where {crfj}jL Q is a set of Jacobi weights. The following result follows from Lemma 3.3 in [2]. 

Lemma 3.1. If v G H m a 4i „(— 1, 1) for some m > 1 and <f> e Vn, then for the Jacobi- Gauss- Lobatto integration, 
we have 

(V,^, - (v,<t>) N \< CN^ || 8™V \\ ua+ m-X,p+ m -l\\ <\> L„,« . 

3.2. Auxiliary results 

By the definitions of the inner product (3.3), the discrete inner product (3.6), and the notations given at the 
beginning of this section, we can rewrite (2.1) at t = t n+ i, (2.5), and (2.6), respectively, as 

*(W)= E^^ + i^^r^iO^iO)),!) a io . (3-7) 



x n +i 

and 

r p - 



[ E^ h) + ^ ) ^)^ IN f^ 1 ) .... ( 3 - 8 ) 



JN 



On the other hand, since each {AI IN,P / n +i(sj)} or { AI IN f n+ i(sj)} is essentially a linear combination of 
parts of {f(ti,Xi)}._ or {f(U, a;j)}™_ and f(t n+ i, x p +1 ), we can also formally rewrite (2.5) and (2.6) as 

l " C *" l ~ 1 t k 1 t " 

Sn+i= J! "l^o* + frr(^y i ) a [X] ai '"+ 1 ^( <i ' a; ^ + a n+i,n+i/(*n+i>a ; n+i) > (3.10) 
fc=0 ' t=0 

and 

<+r = E W + rL ( 2) 8 S^ t1 ' 1 ')' (3 - n) 

where {a^n+i}™^ 1 anc ^ {^i.n+i}"=o are se ^ s °f rea l numbers depending on the number of the interpolating nodes 
IN and the positions of those Jacobi nodes in the interval [0,t n +i]. The formulae (3.10) and (3.11) can help us 
to understand the error analysis we will be performing. First, we have the following estimate. 
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Theorem 3.1. If f(t,x) is sufficiently smooth on a suitable set S € M 2 , and x{t) is also regular enough w.r.t. 
t, then there is a constant C\, independent of n and h, such that 



1 



r(a) Jo 



" (t n+1 - rr- x f{r, x(r))dr - -Lfe)" £ b i<n+1 f{t h x(U)) 



r(a) v 2 



i=0 



a iIN 



< Cit% +1 h 



(3.12) 



Proof. By the definitions of {PI IN ' P f n+1 (s 3 )} and (3.6), we have 



T(a) 



J " +1 (*n+i ~ rf^fir, x(r))dr - -L (*^)° £ b i>n+1 f(U, x z ) 



T(a) v 2 



/„+l (•,£„+!(•)), 1 



-lP7^/ n+1 (.),l 



+ ^/ n+1 (-,x n+1 (.)),lJ ^ ~ ( P/ 7 ^ P /„+i(-),l 

^n + l.l + ^n + 1.2, 



- ( / n+ i(-,£ n+ i (•)),! 



JN. 



JN 



(3.13) 



where 



7„+i,i 

In+1,2 



1 /^n+l\a 



r( a y 2 j 

' / tn+l \ a 



/n+l(-,in+l(-)),l 
/„+! (•,#„+! (•)), 1 



JN 



fn+l(-,X n+ i(-)), 1 
P/™^/ n+1 (.),l 



JJV. 



JJV. 



(3.14) 



Under the assumption of that f(t,x) is sufficiently smooth w.r.t. t, from Lemma 3.1, we know In+i,i can 
be sufficiently small (because m can be arbitrarily large), say, machine accuracy. Using the theories of the 
Lagrange interpolation and of Gauss quadrature [16], we have 



In+l 



1 ,t 



+ 1,2 



n+1 \ a 



L(a) v 2 



JN 



£ w j'( /«+l( s J' i «+l( s j)) _ PlIN ' P fn+l{Sj) 
3=0 



JN 



< C(f,a,IN,JN)t% +1 h IN 



(3.15) 



The inequalities in (3.15) hold because of the positiveness of the Gauss quadratures and J^j^o^j = I-iO- ~ 
s^ds = 2 a /a. 

Next we come to a result corresponding to the corrector formula. Since the proof of the following theorem is 
very similar to the above one, we omit it. 

Theorem 3.2. If f(t,x) is sufficiently smooth on a suitable set S € M 2 , and x{t) is also regular enough w.r.t. 
t, then there is a constant C2, independent of n and h, such that 



r (o0 Jo 



1 ,t 



n+1 



(tn+i - r)"- 1 /^, x(r))dr - frl (^ i ) Q £ a iin+l f(t u x{U)) < C< +1 /* /A 



L(a) 2 



(3.16) 



i=0 
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3.3. Error analysis for the Jacobi-predictor-corrector approach 

In this subsection, we present the main result on the error estimate of the Jacobi-predictor-corrector ap- 
proach, the proof of which is based on the mathematical induction and the results given in the Subsection 
3.2. Through the following result we can observe anther significant advantage of the presented method — the 
convergence order is exactly the number of the interpolating nodes IN, in other words, you can get your desired 
convergent order just by choosing the enough number of interpolating nodes. 



Theorem 3.3. If f(t,x) is sufficiently smooth on a suitable set S E M 2 , h < 1, and x(t) is regular enough 
w.r.t. t, then for the Jacobi-predictor-corrector approach (2.5) and (2.6) and for some suitably chosen T, there 
is a constant C (depends on a, IN and JN), independent of n and h, such that 



max | x(t n+ i) - x n+ i |< Ch , 

l<n+l<W 



(3.17) 



where N = T/h. 



Proof. We use the mathematical induction to prove this theorem. 

a) First we prove Eq. (3.17) holds when n + 1 = IN: Since / is sufficiently smooth, / is legitimately Lipschitz 
continuous. Denoting L as the Lipschitz constant of / w.r.t. its second parameter x, then by (2.1), (3.11) and 
Theorem 3.1, there exists 



| x(t n+1 ) -x n+1 | 
1 



< 



T(a) 
1 



T(a) 



(t n+ i - r) a - x /(r, x{r))d T - (-^) a £ h, n+1 f(U, x 

t n 

(t n+ i - r) a - x /(T, x(r))dr - hn + if(U, x(U)) 

JN , 



3=0 



1 (tn+l\u 



JN 



T(ce) K 2 ; 0<j<JN 
f 



)" max | PI IN > p f n+1 ( Sj )-AI IN ' p f n+1 ( Sj ) | 

( )< j < . / /V ' ■ 



UJj 



3=0 



< fitg^F + "+ 1 C(IN, JN) ■ max | f(t u x(U)) - f{% Xi ) 

1 (Of + 1) 0<i<n 



<- n +a h IN | n b n+l 



C(IN,JN)- max | x(t t ) - Xi \ 

1 [a + Ij Q<i<n 



Cit^h 1 * + C 3 Lt% +1 ■ max | x(t t ) - x t 



0<i<n 



(3.18) 



We have assumed that the starting error maxo<i<„ = /Ar_i | x(ti) — Xi \ is very small (not deteriorating the 
accuracy of the present algorithm), so the first term in the right hand of the last formula in (3.18) plays the 
leading role, thus 

\ x(t n+1 ) - x P +1 \<C 4 t« +1 h IN . (3.19) 
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Combining the above estimate with (2.1), (3.10) and Theorem 3.2, 



I x(t n +l) - Xn+1 I 



1 



< 



T(a) 
1 



(Wi - rr-'fir, x{r))dr - {^T \ £ ^n+if(U, Xi ) + 



a n +l,n+lf(tn+l,Xn+l) 



i=0 
n+1 



(t n+ i - rr-'fir, x{r))dr - {-^) a £ a^ n+1 f (t Zl x(U)) 



i=0 



1 ,t 



n+1 \ a 



v( a y 2 



< C 2 t a n+1 h IN 



JN 



( PI IN f n+1 ( Sj ) - AI IN f n+1 (s 3 

3=0 



1 /tn+l\a 



JN 



T(a) v 2 ; o<j<JN 



max | PI IN f n+1 (sj) ~ AI IN f n+1 ( Sj ) | • ]T 



UJj 



3=0 



< C 2 t* +l h IN + ^," + \, C(IN, JN) 



r(a+l) 

•max{max | f(U,x(t t )) ~ f{t l ,x l ) |, | f (t n +i, x(t n+1 )) - f(t n+ i, x^+i) |} 

< C 2 t% +1 h IN + t'^M lN, JN) ■ max{ max | zft) - |, | x(t n+1 ) - a£ +a |} 

< C 2 t% +1 h IN + C 5 Lt? l+1 • max{ max | x(U) - x t \,C 4 t« +1 h IN } 

0<i<n 

< (C 2 + C 6 Lt" +1 )t" +1 h IN 
— (C2 + CQLt^ N )tf N h IN 

= (C 2 + C 6 i • IN a h a )IN a h a ■ h IN 

< (C 2 + C 6 i • IN a )IN a ■ h IN := C7i 7Ar , 



(3.20) 



where the last inequality holds since h < 1. 

b) We further prove Eq. (3.17) holds for any n+1 > IN: Assume that maxo<i< n +i | x{t n +i)—x n +i \< Ch IN , 
then we are going to prove that | x(t n+2 ) — x n + 2 |< Ch IN . Since the discussions are similar to the ones given 
in a), we briefly present them, 



I x(t n+2 ) - x% +2 | 
< Crt" 2 h IN + C 3 Lt? l+2 • max | x(U) - Xl \ 



0<i<n+l 



< dT a h 1N + C 3 LT a • max I x(t t ) - Xi 

0<i<n+l 



< (Ci + C 3 CL)T a h 



n in 



(3.21) 



and 



x(t n+2 ) - x n+2 

< C 2 t^ +2 h IN + C 5 Lt% +2 • max{ max | x(U) -Xi\, \ x(t n+2 ) - x% +2 } 

U<?<n+1 



< C 2 T a h IN + C 5 LT a • max{C7i Jiv , {d + C 3 CL)T a h 1N }- 



IN 



JNi 



(3.22) 
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by choosing T sufficiently small, we can make sure that C^CLT a ,CiT a ,C§CLT a , as well as C^T" are all 
bounded by C/2, thus 

I x(t n+2 ) - X n+2 | 

< C 2 T a h IN + C 5 LT a ■ Ch IN 

< Ch IN . (3.23) 



Remark 3.1. In practical computations, the Improved-Adams' methods proposed in [4] can be used to start the 
algorithm. We can take the step-length ho discussed in the algorithm description in Section 2 as h - 10~ fc , where 
k > 1 is a given integer, then by the result in [4], there exists 

max I x(ti) — Xi 

0<i<IN-X 



= O(h 



min{l+2a,2} 


10 -fe(i+2a) . o(h 1+2a ), ifO < a < 0.5 



-2* nfu2, «„^ n * ( 3 - 24 ) 



10" 2fe ■ 0(h 2 ), if a > 0.5. 



If taking ft, = 10 m , where m is a given positive integer, by a simple computation, we obtain that \ x(t n+ i) — 
x n +i \— 0(h IN ), as long as the integers k and m satisfy 

IN<{(} + £K + & f <^°- 5; (3.25) 
12 + —, if a > 0.5. v ; 

4. Modifications of the Algorithm 

We have completed the description of the basic algorithm and its most important properties. The convergent 
order of the algorithm is exactly equal to the number of the interpolating points. As is well known, in practice, 
the effectiveness of the algorithm is closely related to the regularity of the solution of the equation. For the 
fractional ordinary differential equation, most of the time, the smoothness of its solution at t = is weaker than 
other places [6] . So we will simply discuss this issue in the following. Another issue we will also mention is how 
to use this algorithm when a is very close to 0. 



4.1. The function / or x is not very smooth at the starting point 



When the smoothness of / or x is weaker at the initial time zero than other time, the sensible way is to divide 
the interval [0,T] into two parts [0, To] and [Tq,T], where To is a small positive real number. For the small 
interval [0, To], we use the Gauss-Lobatto quadrature with the weight function oj(s) = 1. For the remaining 
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part [To,T], the algorithm provided above is employed. That is, 



11 



fe=0 



1 /■*»+ 



t»+i 



E ^ + ^ / (Wi-rr-VMr))^ 



fc=0 



r (") Jo 



-1 
IN 



k \ *o "l»^ 

k=0 y ' j=0 



1 ftn+1 —Tg^a 



J=0 



(4.1) 



where JjV, {Qj}^ and {tj-}^) correspond to the number of, the weights of, and the values of the Gauss- 

Lobatto nodes with the weight uj(s) = 1 in the interval [0, To], respectively. The values of {/(t, , x(rj)) }^_ can 
be computed as in the starting procedure. Since / and x are continuous in the interval [0,To], by the theory of 
Gauss quadrature [16] and the analysis above, we can see that if JN is a big number the accuracy of the total 
error still can be remained. 



4.2. The value of a is very small 

In this subsection, we discuss the case that a is very small, although it happens very less often. When a 
becomes bigger, the weight of the Gauss-Lobatto quadrature at the endpoint of the righthand side becomes 
smaller, the provided algorithm becomes more robust. When a is very small, say, smaller than 0.1, the weight 
at the endpoint of the righthand side of the interval is much bigger than other places (see the Appendix), which 
may impacts the robustness of the algorithm. There are two choices to deal with this problem: one is to try to 
avoid using the high order interpolation in the algorithm; one is to divide the interval [0, T] into two subintervals 
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[0, Ti] and [Ti,T], then similarly do what we do in the last subsection, that is, 



»(*)= E ^o fc) + r^y/ V^/Mr))* 



fc=0 

1 



r(a) 



(t-T^-VCr.xCr))^ 



W_1 t* r« 1 ' Tl 



E <M' 



fe=0 



H*" + !>)/„ («-^V(r.*(r))* 



1 ft-T^ a 



1 



r(a) 



1 /t-Ti 



r(a) 



5>i/(*i,5(»j)), (4.2) 



3=0 



where JiV, {wj}^, and {tj}^ are the same as those defined in the last subsection. And {f(rj, x ( T j)) } J — 
in the second term of the righthand side of the last formula can be computed by interpolation. 

5. Numerical Experiments 

In this section we present two numerical examples to verify the convergent orders derived above and to 
demonstrate the robustness of the provided methods for big T. We only consider the examples with < a < 2, 
since the algorithm will be more robust for a > 2. All numerical computations were done in Matlab 7.5.0 on a 
normal laptop with 1GB of memory. 

5.1. Verification of the error estimates 

First, we use the following example to verify the convergent order: 

D?x(t) = -x(t) + Y^^t 8 -" + 3 ■ Y^^t 7 ~ a +t 8 ~ 3* 7 , (5-1) 

with the initial condition(s) x(0) = (and x'(Q) = if 1 < a < 2). The exact solution of this initial value 
problem is 

x(t) = t 8 + it 7 . (5.2) 
We start the procedure with the Improved- Adams' methods in [4] as discussed in Remark 3.1, i.e. the values 
of x{t) at to,tx, ■ • ■ , tiN-i are computed by the Improved- Adams' methods. The convergent orders are verified 
at T = 1, and the number of the Jacobi nodes is taken as JN + 1 = 27. The number of the interpolating nodes 
IN is respectively taken as 2, 3, 4 and 5 to expect that the corresponding convergent order is also 2, 3, 4 and 
5. The numerical results of the maximum errors for the Jacobi-predictor-corrector approach are showed in the 
following tables, where "CO" means the convergent order. The nodes and weights used in the Gauss-Lobatto 
quadrature w.r.t. the weig ht functions w(s) = (1 - s)"^ 1 + s)° for various a are given in Appendix. Form the 
results in Table 1 to Table 4, we can see that the data confirm the theoretical results provided in Section 3.3. 

Table 5 and Table 6 give the numerical results of the maximum errors for the fractional Adams methods 
in [8] and for the Improved Adams methods in [4]. The compare of Table 1 to Table 4, Table 5 and Table 6 
indicates that the Jacobi-predictor-corrector approach is effective. 
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See the results in Table 3 and Table 4 for a = 0.1, IN = 4 and 5; like what we discussed in Section 4.2, the 
results tell us that we must be more careful to use the provided algorithm when a is very small (letting IN 
be small or dividing the original interval into subintervals) . However, we still confirm the convergent order by 
taking small T (T = 0.1 and T = 0.001) in Table 7. 



Table 1. The maximum errors for (5.1) when t e [0, 1] and IN = 2. 



v, 

11 


OL — U. 1 


PD 


Q — U.O 


pn 


Q — U.O 


pn 


OL — U. ( 


pn 


1/10 


1.34 le+0 


- 


4.21 le-1 


- 


2.27 le-1 


- 


1.68 le-1 


- 


1/20 


3.75 le-1 


1.84 


8.57 le-2 


2.30 


4.65 le-2 


2.29 


3.99 le-2 


2.07 


1/40 


8.40 le-2 


2.16 


1.70 le-2 


2.34 


9.99 le-3 


2.22 


9.14 le-3 


2.13 


1/80 


1.71 le-2 


2.30 


2.98 le-3 


2.51 


1.89 le-3 


2.40 


2.14 le-3 


2.10 


1/160 


3.74 le-3 


2.19 


6.67 le-4 


2.16 


4.17 le-4 


2.18 


5.84 le-4 


1.87 


1/320 


6.97 le-4 


2.42 


9.84 le-5 


2.76 


1.06 le-4 


1.98 


1.28 le-4 


2.19 


1/640 


1.49 le-4 


2.22 


2.75 le-5 


1.84 


2.78 le-5 


1.93 


3.19 le-5 


2.00 


1/1280 


3.67 le-5 


2.03 


7.84 le-6 


1.81 


8.13 le-6 


1.77 


9.03 le-6 


1.82 


1/2560 


9.30 le-6 


1.98 


2.08 le-6 


1.91 


1.94 le-6 


2.07 


2.40 le-6 


1.91 


h 


q = 0.9 


CO 


Q = 1.2 


CO 


Q = 1.5 


CO 


a = 1.8 


CO 


1/10 


1.51 le-1 




1.47 le-1 




1.47 le-1 




1.49 le-1 




1/20 


4.03 le-2 


1.91 


4.16 le-2 


1.83 


3.95 le-2 


1.90 


3.75 le-2 


1.99 


1/40 


8.84 le-3 


2.19 


8.26 le-3 


2.33 


9.01 le-3 


2.13 


1.08 le-2 


1.80 


1/80 


2.32 le-3 


1.93 


2.05 le-3 


2.01 


2.35 le-3 


1.94 


2.58 le-3 


2.06 


1/160 


6.16 le-4 


1.91 


5.07 le-4 


2.02 


6.14 le-4 


1.94 


6.67 le-4 


1.95 


1/320 


1.48 le-4 


2.05 


1.44 le-4 


1.82 


1.64 le-4 


1.91 


1.67 le-4 


1.99 


1/640 


3.84 le-5 


1.95 


3.91 le-5 


1.88 


3.77 le-5 


2.12 


4.21 le-5 


1.99 


1/1280 


9.75 le-6 


1.98 


1.09 le-5 


1.84 


1.01 le-5 


1.90 


1.11 le-5 


1.93 


1/2560 


2.46 le-6 


1.99 


2.69 le-6 


2.02 


2.73 le-6 


1.89 


2.86 le-6 


1.95 



Taking a — 0.5 and h = 1/40, a comparison of the CPU time needed to solve (5.1) for the fractional 
Adams methods in [8], the Improved Adams methods in [4] and the Jacobi-predictor-corrector approach here 
when IN = 2, 3, 4, 5, is reported in Fig. 1. Fig. 1 illustrates that the computational cost of the Jacobi- 
predictor-corrector approach is O(N). Also notice that, as expected, both the fractional Adams methods and 
the Improved Adams methods exhibit 0(N 2 ) computational complexity. 

Although, from Fig. 1, the Jacobi-predictor-corrector approach takes more time to reach the terminate 
time T, when T is small, for example, when T = 1, the CPU time of the fractional Adams methods, the 
Improved Adams methods and the Jacobi-predictor-corrector approach when IN = 2, 3, 4, 5, respectively are 
0.078, 0.109, 0.203, 0.281, 0.313, and 0.359 (sec). While, by Table 1 to Table 6, we can see that, when N = 160, 
the maximum errors of the Jacobi-predictor-corrector methods, 9.99 le-3, 7.19 le-4, 1.43 le-5, 1.08 le-6, are 
much smaller than those of the other two methods, 4.65 le-2, 2.95 le-2. 

Table 8 shows the CPU time (sec) and the steps N needed to solve (5.1) when a — 0.5 with the maximum 
error 1.0 x 10~ 3 , for the fractional Adams methods in [8], the Improved Adams methods in [4] and the Jacobi- 
predictor-corrector approach here when IN — 2, 3, 4, 5. The consumed CPU time presented in Table 8 shows 
that the fractional Adams methods generates the numerical solution with the same accuracy as the other two 
methods, but uses much less CPU time. This advantage is more obvious as the terminate time goes long. It 
further demonstrates the quickness and efficiency of the Jacobi-predictor-corrector method. 
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Table 2. The maximum errors for (5.1) when t G [0, 1] and IN = 3. 



h 


a = 0.1 


CO 


a = 0.3 


CO 


a = 0.5 


CO 


a = 0.7 


CO 


i / 1 a 
1/10 


C A *7 1 n 1 

b.4( le-l 




1 CA 1 n 1 

1.50 le-l 




a ad i n o 
D.D9 le-2 




A 1 n O 

4.2b le-2 




i /on 
1/2U 


O K A 1 n O 

o.o4 Ie-2 


O AO 


1 CA 1 n O 

1.59 le-2 


O O A 

6. 24 


7 A1 1 n O 

(.Ul le-j 


o 

3.2o 


C OO 1 n o 

0.2,3 le-j 


O AO 

3.U3 


i /in 
1/4U 


A £i /I 1 n O 


QIC 

3.15 


1 C7 1 n O 

l.o I le-o 


O O A 

3.34 


7 1 A 1^ /I 

/.ly le-4 


o oo 
3.28 


C 70 1 n A 

o.l I le-4 


O 1 A 

3.19 


i /on 
1/80 


1 AA 1 n O 

I.UU le-3 


O 07 

6.11 


1 97 1 ^ A 

1.3 r le-4 


O CO 

3.02 


a o c in £ 

b.8o le-5 


O OA 

3.39 


7 1 A 1 n c 

(.19 le-o 


O AA 

2.yy 


1 / 1 C. A 

1/ IbU 


1 AT 1 n /I 

l.U( le-4 


o oo 
0.22 


1 OA 1 n F. 

1.39 le-o 


O O A 

3.3U 


7 AC 1 n c 

(.Uo le-o 


o oo 
3.28 


A 77 1 n 

9. ( ( le-o 


o oo 
2.88 


1 /ooa 
I/ozU 


1 AO 1 n C 

1.U2 le-o 


O OA 

3.39 


1 AC 1 n C 

1.06 le-6 


O 71 

3. ( 1 


A CA 1 n 7 

9.oU le-( 


O O A 

2.89 


1 AC 1 n CX 

1.U5 le-o 


o oo 
6.11 


1 /t? /I A 

l/b4U 


i.i ( le-b 


3.13 


1 /I K 1 n 7 

1.45 le-( 


O 07 

2.01 


1 OF 1 n 7 

1.25 le-( 


O AO 

2.yj 


1 Ofi 1 n 7 

1.3b le-( 


OAK 

2.yt> 


1/1280 


1.45 le-7 


3.00 


1 Q>S 1p-8 


2.90 


1.76 le-8 


2.82 




2.78 


1/2560 


1.83 le-8 


2.99 


2.57 le-9 


2.92 


2.17 le-9 


3.02 


2.45 le-9 


3.01 


h 


q = 0.9 


CO 


a = 1.2 


CO 


a = 1.5 


CO 


a = 1.8 


CO 


1/10 


3.51 le-2 




3.27 le-2 




3.24 le-2 




3.27 le-2 




1/20 


4.94 le-3 


2.83 


4.84 le-3 


2.76 


4.46 le-3 


2.86 


4.24 le-3 


2.95 


1/40 


5.40 le-4 


3.19 


5.29 le-4 


3.19 


5.89 le-4 


2.92 


6.79 le-4 


2.64 


1/80 


7.32 le-5 


2.88 


6.70 le-5 


2.98 


7.95 le-5 


2.89 


8.04 le-5 


3.08 


1/160 


9.71 le-6 


2.91 


9.10 le-6 


2.88 


1.05 le-5 


2.92 


1.09 le-5 


2.88 


1/320 


1.23 le-6 


2.99 


1.19 le-6 


2.93 


1.37 le-6 


2.94 


1.38 le-6 


2.98 


1/640 


1.51 le-7 


3.02 


1.65 le-7 


2.86 


1.52 le-7 


3.17 


1.68 le-7 


3.03 


1/1280 


2.08 le-8 


2.86 


2.14 le-8 


2.95 


2.11 le-8 


2.85 


2.18 le-8 


2.95 


1/2560 


2.49 le-9 


3.07 


2.69 le-9 


2.99 


2.78 le-9 


2.92 


3.04 le-9 


2.84 



5.2. Robustness of the algorithm 

Here we study the following equation as an example to show the robustness of the algorithm, 

D?x(t) = -x(t), x(0) = 1, x'(0) = (if 1< a < 2). (5.3) 
It is well known that the exact solution of (5.3) is 

x(t) = E a {-t a ), (5.4) 

where 

00 k 

= E f? fe + iv fle ( a ) >0 ' zeC ' ( 5 - 5 ) 

fc=o ^ a ' ' 

is the Mittag-Leffler function of order a. It is obvious that neither x(t) nor D"x(t) has a bounded first (second) 
derivative at t = when < a < 1 (1 < a < 2), so we deal with (5.3) as we discussed in Section 4.1. Here 
we take Tq = 0.1, JN = 26, JN = 2JN, where Tq, JN and JN are defined as in Section 4.1, and the exact 
solutions are calculated using the function "mlf.m" [19]. The convergent order is also simply verified in Table 
9 and Table 10. 

Further we compute (5.3) with a big time interval, T = 50; the parameters IN, JN, JN are taken the same 
as the above ones and h is taken as 1/10. The exact solutions and relative errors are shown in Fig. 2 with 
a = 0.2 and Fig. 3 with a — 0.5. It can be seen that the relative errors in the interval are less than O(10 -4 ) 
when time is long, which suggests that our method is suitable for the long-time computation. 
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Table 3. The maximum errors for (5.1) when t G [0, 1] and IN = 4. 



h 


a = 0.1 


CO 


a = 0.3 


CO 


a = 0.5 


CO 


a = 0.7 


CO 


i / 1 a 
1/1U 


O OA 1 ~ 1 

2.o9 le-1 




A *71 1 ~ O 

4. / 1 le-2 




1 AO 1 O 

1.4s le-2 




A CA 1 „ O 

4.o0 le-o 




i /on 
1/2U 


i f; i i ~ o 

l.ol le-2 


O AA 

3.9U 


O A 1 1 ~ O 

2.41 le-o 


A OO 

4.28 


C AA 1 A 

0.09 le-4 


A QG 
4.86 


O OO 1 „ c 

8.28 le-o 


o. (b 


i /An 
1/4U 


1 AO 1 ~ O 

l.Uo le-o 


O OA 

3.89 


1 AK 1^/1 

l.Uo le-4 


A CO 

4.02 


1.4o le-o 


Kit; 
0.16 


O *7A 1 ~ 

8. <9 le-o 


O O A 

0.24 


i /on 
1/SO 


O AC 1 A 

o.Oo le-4 


1 oo 

l.so 


/I OO A Ci 

4. 22 le-o 


A G A 
4.D4 


O AA 1 ~ *7 

2.yy ie-< 


0.0 1 


o. m le- / 


O OO 

0.33 


1 / 1 C. A 

1/ 16U 


/I OA 1 ~ O 

4.00 le-o 


O OO 


1 OA 1 ~ ^7 

l.sU le-i 


A EC 
4.00 


1 70 1^,0 

1. (o le-o 


AAA 
4.11 


*7 OA 1 ~ O 

/.2U le-e 


O C1 

0.61 


1 /ooa 
1/ozU 


a o i i^i 

o.ol le-1 


*7 OA 


A CIT 1 ~ A 

4.y< le-y 


C 1 o 

o.ls 


O AO 1 ~ A 

o.92 le-y 


O 1 c 

2.1o 


C 1 A 1 ~ A 

o.lU le-y 


O OO 

o.e2 


1 /t? /I A 

1/64U 


1 A /I 1 „ i 1 

1.94 le+l 


A AK 


O Q7 1 ~ 1 A 

2.sY le-lU 


AAA 
4.11 


O OO 1 „ 1 A 

2.oo le-lU 


A f\A 

4.U4 


O OR A ~ 1 A 

o.3b le-10 


O AO 

0.92 


1/1280 


9 97 1 p+4 


-10.2 


1.57 le-11 


4.19 


1 1e-1 1 

_1_ . C/O 1C A- A- 


3.62 


9 26 1 p-1 1 


3.89 


1/2560 


1.15 le+12 


-25.6 


1.04 le-12 


3.92 


1.03 le-12 


4.23 


1.47 le-12 


3.94 


h 


a = 0.9 


CO 


a = 1.2 


CO 


a = 1.5 


CO 


a = 1.8 


CO 


1/10 


7.94 le-4 




2.82 le-3 




4.34 le-3 




5.00 le-3 




1/20 


2.02 le-4 


1.97 


3.29 le-4 


3.10 


3.65 le-4 


3.57 


3.80 le-4 


3.72 


1/40 


1.59 le-5 


3.67 


2.06 le-5 


4.00 


2.20 le-5 


4.06 


2.69 le-5 


3.82 


1/80 


1.22 le-6 


3.70 


1.31 le-6 


3.97 


1.54 le-6 


3.83 


1.57 le-6 


4.10 


1/160 


9.00 le-8 


3.76 


8.89 le-8 


3.89 


1.08 le-7 


3.84 


1.05 le-7 


3.90 


1/320 


5.74 le-9 


3.97 


5.58 le-9 


3.99 


6.62 le-9 


4.03 


6.78 le-9 


3.95 


1/640 


3.85 le-10 


3.90 


3.94 le-10 


3.82 


3.86 le-10 


4.10 


4.35 le-10 


3.96 


1/1280 


2.36 le-11 


4.03 


2.78 le-11 


3.82 


2.58 le-11 


3.90 


2.70 le-11 


4.01 


1/2560 


1.50 le-12 


3.98 


1.70 le-12 


4.04 


1.67 le-12 


3.95 


1.80 le-12 


3.91 



6. Conclusion 

We provide the Jacobi-predictor-corrector approach for the fractional ordinary differential equations; the 
basic idea is to take the Riemann-Liouville integral kernel as the Jacobi-weight function, and to realize the 
algorithm by doing the Jacobi-Gauss-Lobatto quadrature and polynomial interpolation. The convergent order 
is exactly equal to the number of interpolating nodes IN. The computational complexity is O(N) for a G (0, oo), 
where N is the total computational steps. This is the striking feature/advantage of the algorithm, since the 
computational complexity of numerically solving the fractional ordinary differential equation usually is 0(N 2 ), 
caused by its nonlocal property; when a G (0, 2), it is possible to reduce the computational cost to 0(N log N) 
by combining the short memory principle. 
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Table 4. The maximum errors for (5.1) when t G [0, 1] and IN = 5. 



h 


a = 0.1 


CO 


a = 0.3 


CO 


a = 0.5 


CO 


a = 0.7 


CO 


i /in 
1/10 


A *7C\ 1 ~ 1 

4. 1 y le-2 




i Jn 1. n 
1.49 le-2 




A n*7 1 ~ o 

4.U/ le-3 




1 OO 1 „ o 

1.2o le-o 




i /on 
1/2U 


l.ol le-2 


1.0/ 


O C /I 1 ~ /I 

3.54 le-4 


c a n 
0.4U 


7 OO 1 C 

i\22 le-o 


C OO 

0.82 


i ^ n i c 
1.4y le-o 


0.3 1 


i /in 
1/4U 


Oil 1 ~ 1 

a. 11 le-l 


-5. to 


i .14: le-6 


C CO 

0.52 


i no i a 

l.Us le-o 


p. n*7 


O *7 A 1 „ *7 

2. /4 le-< 


o. (b 


i /on 
1/SO 


1 O C 1 ~ /I 

l.zo le+4 


ion 
-13.9 


1 CO 1^7 


C Q~\ 

O.Ol 


i on i ~ o 
l.o9 le-o 


CX OO 

6.28 


1 RQ 1 „ O 

l.oo le-8 


a no 
4.U3 


1 / 1 C. A 

1/ 160 






001 i ~ n 

3.31 le-9 


ceo 
0.08 


i no i ~ 1 n 
l.yo le-lU 


6.1 / 


*7 at 1 ~ in 

r.4f le-10 


A ACi 

4.4y 


1/320 








5.99 




3.31 


9 57 1 e-1 1 


4.86 


1/640 






1.60 le-12 


5.02 


5.64 le-13 


5.11 


8.92 le-13 


4.85 


1/1280 






5.06 le-14 


4.98 


1.87 le-14 


4.92 


3.06 le-14 


4.86 


h 


a = 0.9 


CO 


a = 1.2 


CO 


a = 1.5 


CO 


a = 1.8 


CO 


1/10 


2.69 le-4 




4.47 le-4 




7.36 le-4 




8.87 le-4 




1/20 


1.49 le-5 


4.18 


2.60 le-5 


4.10 


2.94 le-5 


4.64 


3.14 le-5 


4.82 


1/40 


6.52 le-7 


4.51 


8.53 le-7 


4.93 


9.18 le-7 


5.00 


1.10 le-6 


4.84 


1/80 


2.47 le-8 


4.72 


2.63 le-8 


5.02 


3.25 le-8 


4.82 


3.07 le-8 


5.16 


1/160 


9.01 le-10 


4.78 


9.88 le-10 


4.73 


1.14 le-9 


4.84 


1.08 le-9 


4.83 


1/320 


2.90 le-11 


4.96 


2.86 le-11 


5.11 


3.31 le-11 


5.10 


3.42 le-11 


4.98 


1/640 


9.53 le-13 


4.93 


1.01 le-12 


4.83 


1.00 le-12 


5.05 


1.10 le-12 


4.96 


1/1280 


3.46 le-14 


4.78 


3.38 le-14 


4.90 


3.24 le-14 


4.95 


3.46 le-14 


4.99 



Table 5. The maximum errors of fractional Adams methods for (5.1) when t G [0, 1 



h 


Q = 0.1 


CO 


a = 0.3 


CO 


a = 0.5 


CO 


a = 0.7 


CO 


1/10 


2.09 le+0 




8.45 le-l 




4.51 le-l 




2.89 le-l 




1/20 


1.17 le+0 


0.83 


3.32 le-l 


1.35 


1.46 le-l 


1.63 


8.14 le-2 


1.83 


1/40 


5.64 le-l 


1.06 


1.23 le-l 


1.43 


4.65 le-2 


1.65 


2.28 le-2 


1.83 


1/80 


2.49 le-l 


1.18 


4.53 le-2 


1.44 


1.50 le-2 


1.64 


6.47 le-3 


1.82 


1/160 


1.06 le-l 


1.23 


1.68 le-2 


1.43 


4.90 le-3 


1.61 


1.85 le-3 


1.80 


1/320 


4.52 le-2 


1.24 


6.30 le-3 


1.41 


1.63 le-3 


1.59 


5.38 le-4 


1.79 


1/640 


1.92 le-2 


1.23 


2.39 le-3 


1.40 


5.51 le-4 


1.57 


1.58 le-4 


1.77 


h 


a = 0.9 


CO 


a = 1.2 


CO 


a = 1.5 


CO 


a = 1.8 


CO 


1/10 


2.16 le-l 




1.72 le-l 




1.57 le-l 




1.53 le-l 




1/20 


5.55 le-2 


1.96 


4.17 le-2 


2.04 


3.81 le-2 


2.04 


3.75 le-2 


2.03 


1/40 


1.42 le-2 


1.97 


1.02 le-2 


2.04 


9.35 le-3 


2.03 


9.28 le-3 


2.01 


1/80 


3.65 le-3 


1.96 


2.49 le-3 


2.03 


2.31 le-3 


2.02 


2.31 le-3 


2.01 


1/160 


9.39 le-4 


1.96 


6.11 le-4 


2.03 


5.71 le-4 


2.01 


5.75 le-4 


2.00 


1/320 


2.42 le-4 


1.95 


1.50 le-4 


2.02 


1.42 le-4 


2.01 


1.44 le-4 


2.00 


1/640 


6.26 le-5 


1.95 


3.70 le-5 


2.01 


3.53 le-5 


2.01 


3.59 le-5 


2.00 
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Table 6. The maximum errors of Improved-Adams methods for (5.1) when t 6 [0, 1]. 
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h 


a = 0.1 


CO 


a = 0.3 


CO 


a = 0.5 


CO 


a = 0.7 


CO 


1/10 


2.05 le+0 




7.68 le-1 




3.74 le-1 




2.29 le-1 




1/20 


1.13 le+0 


0.87 


2.71 le-1 


1.50 


1.01 le-1 


1.89 


5.32 le-2 


2.11 


1/40 


5.21 le-1 


1.11 


8.84 le-2 


1.61 


2.59 le-2 


1.96 


1.21 le-2 


2.14 


1/80 


2.21 le-1 


1.24 


2.81 le-2 


1.65 


6.52 le-3 


1.99 


2.76 le-3 


2.13 


1/160 


9.04 le-2 


1.29 


8.90 le-3 


1.66 


1.63 le-3 


2.00 


6.37 le-4 


2.12 


1/320 


3.68 le-2 


1.30 


2.82 le-3 


1.66 


4.05 le-4 


2.01 


1.49 le-4 


2.10 


1/640 


1.50 le-2 


1.30 


9.00 le-4 


1.65 


1.01 le-4 


2.01 


3.52 le-5 


2.08 


h 


q = 0.9 


CO 


a = 1.2 


CO 


a = 1.5 


CO 


a = 1.8 


CO 


1/10 


1.74 le-1 




1.50 le-1 




1.48 le-1 




1.49 le-1 




1/20 


3.93 le-2 


2.14 


3.56 le-2 


2.08 


3.61 le-2 


2.03 


3.68 le-2 


2.02 


1/40 


9.08 le-3 


2.11 


8.69 le-3 


2.04 


8.96 le-3 


2.01 


9.18 le-3 


2.00 


1/80 


2.15 le-3 


2.08 


2.15 le-3 


2.02 


2.24 le-3 


2.00 


2.29 le-3 


2.00 


1/160 


5.20 le-4 


2.05 


5.35 le-4 


2.01 


5.58 le-4 


2.00 


5.73 le-4 


2.00 


1/320 


1.27 le-4 


2.03 


1.34 le-4 


2.00 


1.40 le-4 


2.00 


1.43 le-4 


2.00 


1/640 


3.15 le-5 


2.02 


3.34 le-5 


2.00 


3.49 le-5 


2.00 


3.58 le-5 


2.00 



Table 7. The maximum errors for (5.1) when a = 0.1, IN = 4 and 5. 



fl 


T = 0.1, IN = 4 


CO 


T = 0.001, IN = 5 


CO 


T/10 


9.64 le-9 




1.45 le-23 




T/20 


6.61 le-10 


3.87 


4.97 le-25 


4.87 


T/40 


3.91 le-11 


4.08 


1.46 le-26 


5.09 


T/80 


2.14 le-12 


4.19 


4.01 le-28 


5.18 


T/160 


1.16 le-13 


4.21 


1.11 le-29 


5.18 


T/320 


5.83 le-15 


4.32 


2.88 le-31 


5.26 


T/640 


3.19 le-16 


4.19 


8.19 le-33 


5.14 



Table 8. The CPU time (sec) and the steps N needed to solve (5.1) when a = 0.5 with the 
maximum error 1.0 x 10 , for the fractional Adams methods in [8], the Improved Adams 
methods in [4] and the Jacobi-predictor-corrector approach here when IN — 2, 3, 4, 5. 



terminal time 



methods 


T = 0.5 


T = 1.0 




T = 1.5 




T = 2.0 






N 


CPU time (sec) 


N 


CPU time (sec) 


N 


CPU time (sec) 


N 


CPU time (sec) 


fractional Adams 


17 


6.25 le -2 


432 


5.72 le+0 


3240 


3.20 le+2 


14200 


6.26 le+3 


Improved Adams 


14 


3.13 le -2 


204 


7.97 le -1 


945 


1.41 le+1 


2831 


1.29 le+2 


Jacobi- 


IN — 2 


11 


9.34 le -2 


119 


6.09 le -1 


492 


2.52 le+0 


1456 


7.75 le+0 


predictor- 


IN — 3 


7 


1.56 le -2 


34 


1.41 le -1 


89 


5.16 le -1 


117 


1.13 le+0 


corrector 


IN — 4 


5 


1.56 le -2 


18 


7.81 le -2 


34 


1.72 le -1 


51 


2.97 le -1 


methods 


IN = 5 


5 


3.13 le -2 


13 


4.69 le -2 


23 


1.09 le -1 


33 


1.88 le -1 
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Jacobi-predictor-corrector, IN=2 




Jacobi-predictor-corrector, IN=3 




Jacobi-predictor-corrector, IN=4 




Jacobi-predictor-corrector, IN=5 




improved Adams 




fractional Adams 









5 10 15 20 



FIGURE 1. The CPU time needed to solve (5.1) when a = 0.5, h = 1/40 for the fractional 
Adams methods in [8] , the Improved Adams methods in [4] and the Jacobi-predictor-corrector 
approach here when IN = 2, 3, 4, 5. 

Table 9. The maximum errors for (5.3) when t £ [0, 1.1] and IN = 2. 



h 


a = 0.2 


CO 


a = 0.5 


CO 


a = 1.2 


CO 


a = 1.8 


CO 


1/10 


4.84 le-3 




2.30 le-3 




1.20 le-4 




3.72 le-4 




1/20 


1.49 le-3 


1.70 


5.10 le-4 


2.17 


3.47 le-5 


1.79 


1.06 le-4 


1.81 


1/40 


4.04 le-4 


1.88 


1.02 le-4 


2.33 


7.83 le-6 


2.15 


2.62 le-5 


2.02 


1/80 


9.96 le-5 


2.02 


1.89 le-5 


2.43 


1.87 le-6 


2.07 


6.35 le-6 


2.05 


1/160 


2.44 le-5 


2.03 


3.95 le-6 


2.26 


5.41 le-7 


1.79 


1.62 le-6 


1.97 



Table 10. The maximum errors for (5.3) when t E [0, 1.1] and IN = 3. 



h 


a = 0.2 


CO 


a = 0.5 


CO 


a = 1.2 


CO 


a = 1.8 


CO 


1/10 


2.77 le-3 




7.30 le-4 




1.11 le-5 




1.64 le-5 




1/20 


6.04 le-4 


2.20 


1.14 le-4 


2.68 


3.23 le-6 


1.78 


3.00 le-6 


2.45 


1/40 


1.06 le-4 


2.51 


1.43 le-5 


2.99 


5.48 le-7 


2.56 


4.64 le-7 


2.69 


1/80 


1.29 le-5 


3.04 


1.40 le-6 


3.36 


8.88 le-8 


2.63 


5.91 le-8 


2.97 


1/160 


1.36 le-6 


3.25 


3.78 le-8 


5.21 


1.09 le-8 


3.03 


7.84 le-9 


2.92 
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FIGURE 2. (a) Exact solution of (5.3) for a = 0.2; (b) Relative errors with a = 0.2, h = 1/10, 
IN = 2 or 3. 





(a) 



(b) 



FIGURE 3. (a) Exact solution of (5.3) for a — 0.5; (b) Relative errors with a = 0.5, h = 1/10, 
IN = 2 or 3. 
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Appendix 

We display the Jacobi-Gauss-Lobatoo nodes and weighs in the reference interval [—1,1] used in the numerical 
experiments in the following tables, where the weight function is (1 — s) a , a = 0.1, 0.3, 0.5, 0.7, 0.9, 1.2, 1.5, 
1.8, respectively, and the number of the quadrature nodes JN + 1 = 27. 

Table 11. 



a = 0.1 



a = 0.3 



nodes 




weights 


nodes 


weights 
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i nnnnnnnnnnnnnnnn 

- 1 .uuuuuuuuuuuuuuuu 


n nm snn/i/i^i 7roi qi 

u.uuiouu^i^iui i oy iy i 


-u.yoyzuiuuzyu4:oyuu 


n 
u 


nnQ748fiQfi8 t; ;i ^111 


n QRQ98*3fi99nzLQ c ;i ^4 
■u.yoyioouiiu'iyji^ 


n niiimii *39/iooi aa 

u.uiiiuii ioz4yyi44 
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n QfizL99fizL*31 71 89^87 


n n9nnm RS9R1 non94 


n Q9/i7n/i8fii ^OQQ^/ifi 
-u.yz^ i u^ooiouyyo^o 


n 
u 
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UZOOOOUUOO^i ( ooio 


H 09 c i97ni 78^81/1(^9 


n 09801 98*3*308^9^9fi 


-U.o / ZUZ { OooOoZ i DUO 


n 
u 


n^fi/l 9*30/1 8*3/1/1890. 


879070^*31 fi*3Q9 c i' : 17 

-u.o I zy ( yooiuoyzoo / 


n n*378/ifi^8^i i 880*3*3 

U.Uo t O^DOoOl looUOO 


-u.ouuuo / yiuiZiZoZiOO 


n 


041 QnQ78zL*38'39n/lfi 
u^tiyuy i o'iooOiU^u 


n snsi nsszi^RQi i sn9 


n nzlfiSI 9S99^!Q7fiS1 n 


79Qfi^M ^4780^79 

-<J. 1 ZiyUOO 1004: t OkjO l A 


n 




n 7*3,1 W I i49 c i48nQ' : iQ^ 


n c i c i89'3099'31 0^*378 


-0.6419886593191621 





0592179803160306 


-0.6445363551226361 


0.0648912240850699 


-0.5450216478217176 





0683915067416097 


-0.5481926424422340 


0.0740349165644332 


-0.4401427133803333 





0780200133502234 


-0.4439511569824626 


0.0832761166251512 


-0.3288753760288977 





0882063337718825 


-0.3333146137264882 


0.0926427678850384 


-0.2128359531749780 





0990774042058216 


-0.2178779134712215 


0.1021706112368920 


-0.0937100816930866 





1107929333210531 


-0.0993051526697154 


0.1119057485004456 


0.0267717676524930 





1235579229632447 


0.0206943649229667 


0.1219082421622975 


0.1468594274088615 





1376412507907670 


0.1403907686322420 


0.1322573016048666 


0.2648084569648291 





1534041057504351 


0.2580585579193416 


0.1430589703422321 


0.3789054830602389 





1713450510950583 


0.3720014765729532 


0.1544578938601842 


0.4874930892080245 





1921744216247991 


0.4805769654678000 


0.1666560217001116 


0.5889938926263092 





2169432942460124 


0.5822198414675918 


0.1799436799255219 


0.6819334593243178 





2472807617819305 


0.6754648613490701 


0.1947540530811646 


0.7649617255263299 





2858641340356679 


0.7589678460750942 


0.2117653286898140 


0.8368726174182031 





3374445171968009 


0.8315250626755174 


0.2321093791361104 


0.8966215953237894 





4113924209918509 


0.8920905904778630 


0.2578499674265770 


0.9433409167504676 





5293048192507985 


0.9397914478954834 


0.2932722023123962 


0.9763527097958942 





7555606909447271 


0.9739404234643343 


0.3493613894433857 


0.9951835446365114 


1 


4176717079984273 


0.9940482649436114 


0.4687789976411144 


1.0000000000000000 


5 


0539800138885518 


1.0000000000000000 


0.4664213940659055 
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a = 0.5 



a = 0.7 



nodes 




weights 


nodes 




weights 


-1.0000000000000000 





0020525595970582 


-1.0000000000000000 





0023401106204443 


-0.9893643555933919 





0126420869565904 


-0.9894438812775599 





0143980191213123 


-0.9644947993303408 





0227208651492852 


-0.9647591646875335 





0258125839599000 


-0.9258270443312571 





0327129904967510 


-0.9263756474471649 





0370189772483234 


-0.8739173433922172 





0425872817104095 


-0.8748413378464713 





0479340538950119 


-0.8095088697871018 





0523089398061103 


-0.8108884559986644 





0584716242404765 


-0.7335232176624819 





0618432836433182 


-0.7354251541586265 





0685470018011806 


-0.6470475078824108 





0711562189729229 


-0.6495229109880948 





0780777884066716 


-0.5513189009359872 





0802144213662633 


-0.5544013837390309 





0869842827280500 


-0.4477069175659654 





0889854705911301 


-0.4514111111711064 





0951897677589023 


-0.3376938533233919 





0974379713316180 


-0.3420143458160455 





1026206972731357 


-0.2228535755897699 





1055416672643979 


-0.2277642960175138 





1092067627458629 


-0.1048290089639729 





1132675500578416 


-0.1102830750128124 





1148808056093876 


o ni /ifioi ^^707^71 
u.ui^oy iod / y / o / ido 


u 


1 90^870^3^9399/10 


0087^1 3980^738^0 


u 


1 1 0^78^91 81 1 7^37 


0.1339976779295242 





1274767027659907 


0.1276787341676305 





1232378787845395 


0.2513831064264158 





1339091080692547 


0.2447807621576336 





1257981209039414 


0.3651683196413840 





1398621532108058 


0.3584048089141069 





1271981640581311 


0.4737254891175540 





1453145279145870 


0.4669376500627146 





1273740449364915 


0.5755015796722314 





1502467141498365 


0.5688383448113328 





1262548372106364 


0.6690405673158345 





1546410560090228 


0.6626601132252578 





1237559426680637 


0.7530042693246988 





1584818229166975 


0.7470708756553147 





1197675894104367 


0.8261914884659855 





1617552659442731 


0.8208721610920136 





1141338609152897 


0.8875551974923557 





1644496670298319 


0.8830161106263974 





1066110258518362 


0.9362175180611009 





1665553809272196 


0.9326203125145781 





0967739548306226 


0.9714822797862176 





1680648697345029 


0.9689801330973330 





0837631042002331 


0.9928449797512120 





1689727298783471 


0.9915771271381129 





0653397616047768 


1.0000000000000000 





0846378557289007 


1.0000000000000000 





0196518498509760 
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a = 0.9 



a = 1.2 



nodes 




weights 


nodes 


weights 


-1.0000000000000000 





0026680954862362 


-1.0000000000000000 


0.0032485813206930 


-0.9895222260184334 





0163990211169060 


-0.9896375857934652 


0.0199364618399667 


-0.9650196167842330 





0293282050572500 


-0.9654031458102729 


0.0355265819180070 


-0.9269161710244489 





0418988082171859 


-0.9277121945697792 


0.0504662938531113 


-0.8757518197224115 





0539656088052377 


-0.8770928511031259 


0.0644948884276532 


-0.8122480503073657 





0653836420874911 


-0.8142509047810206 


0.0773638096423790 


-0.7372998407745818 





0760144717348032 


-0.7400620590862488 


0.0888474771253680 


-0.6519633346267486 





0857281246221285 


-0.6555600135812010 


0.0987484771441300 


-0.5574410233343691 





0944045439109477 


-0.5619221255221779 


0.1069017533569795 


-0.4550648215271042 





1019348490180515 


-0.4604530253538590 


0.1131780974268036 


-0.3462773061846129 





1082224273322169 


-0.3525664467804354 


0.1174869327523327 


-0.2326113924086162 





1131838455110155 


-0.2397655327004557 


0.1197783563983938 


-0.1156687346912121 





1167495593531492 


-0.1236218938939909 


0.1200444135582023 


n 0090098/11070^09 

u.uuzyuzo^iu / udouz 


o 
u 


1 1 88fi/l 30^0701 940 


00^7^*371 *31 91 A 0fi7 


1 1 8^1 Q^Qfififi^^l Q 

u. i looiyoyDDDOOoiy 


0.1214325564431927 





1194877759976240 


0.1121967999388062 


0.1146805830124639 


0.2382502228570740 





1185936294194118 


0.2285862886357495 


0.1092452502518708 


0.3517097756999029 





1161699441834089 


0.3417931453037846 


0.1021710400770435 


0.4602124683760809 





1122178437862630 


0.4502401041519408 


0.0936527803497127 


0.5622293993378423 





1067500287515353 


0.5524162161956613 


0.0839201325248429 


0.6563230543010921 





0997882844987957 


0.6468978998418504 


0.0732349201079911 


0.7411675592086939 





0913594924131627 


0.7323687733098204 


0.0618887495635877 


0.8155673560609422 





0814889890869870 


0.8076379911244389 


0.0502016389583398 


0.8784740302618931 





0701886634133376 


0.8716568246106763 


0.0385230330239028 


0.9290010193400755 





0574330648431000 


0.9235332376746080 


0.0272382282726720 


0.9664358148554766 





0431024991510826 


0.9625441648076186 


0.0167880196253709 


0.9902477964312976 





0268013946767667 


0.9881446058894564 


0.0077265395634170 


1.0000000000000000 





0052794393153508 


1.0000000000000000 


0.0008846215676251 
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a = 1.5 



a = 1. 



nodes 




weights 


nodes 




weights 


-1.0000000000000000 





0039558421325122 


-1.0000000000000000 





0048176566867522 


-0.9897504320058830 





0242407262612330 


-0.9898608459541634 





0294787408291468 


-0.9657783446854108 





0430450442150069 


-0.9661454822450298 





0521665283985357 


-0.9284910158577785 





0608075612568833 


-0.9292531881022244 





0732935305513724 


-0.8784051049048023 





0771197426976410 


-0.8796895020922768 





0922636940579904 


-0.8162111679198959 





0916086237361090 


-0.8181301941668921 





1085561087536052 


-0.7427661984333459 





1039548961583067 


-0.7454140912397981 





1217534509743617 


-0.6590821003162865 





1139030192180424 


-0.6625319256813679 





1315582543475505 


-0.5663118098346140 





1212688126551674 


-0.5706128999526579 





1378033510671248 


-0.4657334305951538 





1259447610428847 


-0.4709093210906420 





1404566089311143 


-0.3587326320143288 





1279029357867882 


-0.3647795460854125 





1396198574132639 


-0.2467835618430933 





1271954603677246 


-0.2536694784396956 





1355220825893431 


-0.1314285381791132 





1239525097131300 


-0.1390928704191656 





1285072066852624 


n m /i^fisni fii 701 ^o 


n 
u 


1 1 8*3770081 Q9fi1 fiS 


n 099(^1 07001 8/1 £10/18 


o 
u 


1 1 QOI fiQQQI 8099^ 

±±yu±oyyy±ouzzou 


0.1031173793924301 





1107424642216554 


0.0941900949493554 





1075698748132034 


0.2190769482282871 





1013752491912932 


0.2097182393627444 





0947365082759497 


0.3320243383882507 





0906530919063302 


0.3223997996747073 





0811133270750175 


0.4404034835594873 





0789886148348753 


0.4306996300173916 





0672950260065822 


0.5427212564354959 





0668171834999999 


0.5331422907681873 





0538472733927563 


0.6375680414418941 





0545831738546201 


0.6283321577931673 





0412807461037388 


0.7236371592440043 





0427259833787319 


0.7149724532569940 





0300275324422679 


0.7997428790488688 





0316662190978989 


0.7918829524158366 





0204207704286011 


0.8648367821754457 





0217924888159768 


0.8580161681069255 





0126781227682169 


0.9180222971008329 





0134491973972607 


0.9124719559264060 





0068892709351844 


0.9585674487321285 





0069256857723223 


0.9545111721119047 





0030068635217618 


0.9859178863652559 





0024466760492113 


0.9835752524825214 





0008386286076415 


1.0000000000000000 





0001742117099051 


1.0000000000000000 





0000387924881512 
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